// ==================================================================================
// Shared Genomics Project MPI Codebase
// Version 1.0 30/04/2010
//
// (c) 2010 University of Manchester all rights reserved
//
// This file is distributed under the GNU General Public License, Version 2.  
// Please see the file COPYING.txt for more details
// ==================================================================================

/*!
\file
\brief	Implementation of P-value Adjustment for multiple testing. 
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "adjusted_2d.h"
#include "copyfile.h"
#include "pmaths.h"
#include "ppair.h"
#include "utils.h"
#include "suffixes.h"

struct adjusted_2d* adjusted_2d_create_array(int n) {
	struct adjusted_2d* values = NULL;

	if (n > 0) {
		values = (struct adjusted_2d*) malloc(sizeof(struct adjusted_2d) * n);
		if (values != NULL) {
			int i;
			for (i = 0; i < n; i++) {
				values[i].bonf  =  values[i].fdr_bh = values[i].holm = 1;
				values[i].coeff = values[i].z = values[i].p = -1;
				values[i].p1 = values[i].p2 = -1;
			}
		}
	}

	return values;
}

int adjusted_2d_calc(struct adjusted_2d *a, int nA) {	
	int i, nS = 0; // nS is elements used in the sp & schi arrays.
	double t;
	BOOL jobDone = FALSE;
	struct ppair *sp = NULL;
	VECTOR_T pv_BH = NULL, pv_holm = NULL;

	if (a == NULL || nA <= 0) return 0;

	// If have storage, do the sums.
	sp = ppair_create_array(nA);

	if (sp != NULL) {
		// Calculate the pvalues.
		for (i = 0; i < nA; i++) {
			if (a[i].p >= 0) {
				struct ppair *pt = &sp[nS++];

				pt->p = a[i].p;
				pt->l = i; 
			}
		}
		qsort(sp, nS, sizeof(struct ppair), ppair_cmp); 

		t = (double) nS;
		pv_BH = create_vector(nA, 1.0);
		pv_holm = create_vector(nA, 1.0);

		// Works only if not zero valid tests computed.
		if (nS > 0 && pv_BH != NULL && pv_holm != NULL) {
			qsort(sp, nS, sizeof(struct ppair), ppair_cmp); 

			// Holm 
			pv_holm[0] = sp[0].p*t > 1 ? 1 : sp[0].p*t;
			for (i = 1; i < t; i++) {
				double x = (nS-i) * sp[i].p < 1 ? (nS - i) * sp[i].p : 1;
				pv_holm[i] = pv_holm[i-1] > x ? pv_holm[i-1] : x;
			}

			// BH
			pv_BH[nS-1] = sp[nS-1].p;
			for (i= nS - 2; i>=0; i--) {
				double x = (t/(double)(i+1)) * sp[i].p < 1 ? (t/(double)(i+1)) * sp[i].p : 1 ;
				pv_BH[i] = pv_BH[i+1] < x ? pv_BH[i+1] : x;
			}

			// Map the adjusted_2d pvalues back onto the source array.
			for (i = 0; i < nS; i++) {
				struct adjusted_2d *pv = &a[sp[i].l];

				pv->bonf = sp[i].p*t > 1 ? 1 : sp[i].p*t;
				pv->holm = pv_holm[i];
				pv->fdr_bh = pv_BH[i];
			}

			jobDone = TRUE;
		}
	} 

	if (sp != NULL) free(sp);
	if (pv_BH != NULL) free(pv_BH);
	if (pv_holm != NULL) free(pv_holm);

	if (jobDone) return 1;
	else return 0;
}

void adjusted_2d_echo(struct adjusted_2d *a) {
	if (a != NULL) fprintf(stdout, "p1=%i, p2=%i, z=%f, p=%f, bonf=%f holm=%f fdr_bh=%f\n", a->p1, a->p2, a->z, a->p, a->bonf, a->holm, a->fdr_bh);
}

int adjusted_2d_write_and_copy(struct adjusted_2d *a, int nA, BOOL *items_to_keep, struct selected_options *ops, int rank) {
	FILE *out = NULL;
	char localFilepath[FILE_PATH_LENGTH], remoteFilepath[FILE_PATH_LENGTH];
	BOOL jobDone = FALSE;
	int nToKeep = 0, nWrote = 0, i;

	if (a == NULL || nA <= 0 || items_to_keep == 0 || ops == NULL || rank <  0) return 0;
	else if (strlen(ops->ldir) == 0 || strlen(ops->rdir) == 0 || strlen(ops->szRunId) == 0) return 0;

	for (i = 0; i < nA ; i++) if (items_to_keep[i]) nToKeep++;
	if (nToKeep == 0) return 0; // Nothing to write to file.

	sprintf(localFilepath, "%s%s%s_%i", ops->ldir, ops->szRunId, SUFFIX_PADJ_FILE, rank);
	sprintf(remoteFilepath, "%s%s%s", ops->rdir, ops->szRunId, SUFFIX_PADJ_FILE);

	out = fopen(localFilepath, "wb");
	if (out != NULL) {
		if (fwrite(&nToKeep, sizeof(int), 1, out) == 1) {
			for (i = 0; i < nA; i++) {
				if (items_to_keep[i]) {
					if (fwrite(&a[i], sizeof(struct adjusted_2d), 1, out)) nWrote++;
					else break;
				}
			}
		}
		fclose(out);
		if (nToKeep == nWrote) {
			if (copyAndDeleteFile(localFilepath, remoteFilepath)) jobDone = TRUE;
		} else 
			deleteFile(localFilepath);
	}
	if (jobDone) return 1;
	else return 0;
}

// QSORT comparator function.
int adjusted_2d_cmp(const void* p1, const void* p2) {
	const struct adjusted_2d *a = (const struct adjusted_2d*) p1;
	const struct adjusted_2d *b = (const struct adjusted_2d*) p2;

	if (a == NULL || b == NULL) return 0;

	if (a->p1 < b->p1) return -1;
	else if (a->p1 > b->p1) return 1;
	else if (a->p1 == b->p1) {
		if (a->p2 < b->p2) return -1;
		else if (a->p2 > b->p2) return 1;
		else if (a->p2 == b->p2) return 0;
	}

	return  0;
}

struct adjusted_2d* adjusted_2d_load(struct selected_options *ops, int *size) {
	FILE *in = NULL;
	struct adjusted_2d* data = NULL;
	char remoteFilepath[FILE_PATH_LENGTH];


	if (ops == NULL || size == NULL) return NULL;
	else if (strlen(ops->rdir) == 0 || strlen(ops->szRunId) == 0) return NULL;

	*size = 0;

	sprintf(remoteFilepath, "%s%s%s", ops->rdir, ops->szRunId, SUFFIX_PADJ_FILE);
	in = fopen(remoteFilepath, "rb");
	if (in != NULL) {
		if(fread(size, sizeof(int), 1, in) == 1) {
			if (size > 0) {
				data = (struct adjusted_2d*) malloc(sizeof(struct adjusted_2d) * (*size));
				if (data != NULL) {
					if (fread(data, sizeof(struct adjusted_2d), *size, in) != *size) {
						free(data);
						data = NULL;
					} else qsort(data, *size, sizeof(struct adjusted_2d), adjusted_2d_cmp);
				}
			}
		}
		fclose(in);
	}

	return data;
}

int adjusted_2d_search(struct adjusted_2d *a, int nA, struct adjusted_2d *x)
{
	int result = 0;
	if (a != NULL && nA > 0) {
		int l, r, c;
		for (l = 0, r = nA-1; l <= r;) {
			c = l + (r - l) / 2;
			result = adjusted_2d_cmp(&a[c], x);
			if (result == 1) r = c - 1; // greater than
			else if (result == -1) l = c + 1; // less than
			else return c;
		}
	}

	return -1;
}